Exact Results for Three-Body Correlations in a Degenerate One-Dimensional Bose Gas 
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Motivated by recent experiments we derive an exact expression for the correlation function enter- 
ing the three-body recombination rate for a one-dimensional gas of interacting bosons. The answer, 
given in terms of two thermodynamic parameters of the Lieb-Liniger model, is valid for all values 
of the dimensionless coupling 7 and contains the previously known results for the Bogoliubov and 
Tonks-Girardeau regimes as limiting cases. We also investigate finite-size effects by calculating the 
correlation function for small systems of 3, 4, 5 and 6 particles. 



Quantum fluctuations are well known to have a pro- 
found influence on the physics of one-dimensional (ID) 
systems. Due to quantum fluctuations, a ID gas of 
weakly interacting bosons does not exhibit true long- 
range order and hence does not undergo Bose-Einstein 
condensation at any temperature. With increasing re- 
pulsion between the bosons, the gas exhibits a smooth 
crossover to the strong coupling regime, where the cor- 
relations in the positions of the particles endow it with 
fermion-like thermodynamic properties. A strong repul- 
sion between the bosons also has a marked effect on the 
local fluctuation properties, leading, for example, to a 
fermion-like suppression of the fluctuations of particle 
density. In this paper we shall focus on a three-body 
correlation function, which is directly related to the life- 
time of a Bose gas unstable with respect to three-body 
recombination processes. It has been studied in a recent 
experiment, where measurements of the three-body re- 
combination rate for trapped atoms were performed Q]. 

The idea of using three-body recombination for a mea- 
surement of local correlations was originally put forward 
by Kagan et al. Q. These authors showed that the low- 
temperature recombination rate for a 3D Bose gas is pro- 
portional to the local three-body correlation function g 3 
defined by the ground-state expectation value 
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Here h(r) = ^'(r)^f(r) is the particle density opera- 
tor, iff' ( 1 I r ) denote boson creation (annihilation) oper- 
ators and the symbol :: stands for normal ordering. It 
was noted that, due to the factor-of-six suppression of 
gz in a Bose-Einstein condensate at zero temperature as 
compared to the non-condensed Bose gas, three-body re- 
combination can be used as a diagnostic tool for distin- 
guishing between the condensed and the non-condensed 
phases. This was later confirmed experimentally 

In the experiment of Ref. Q the recombination process 
was used for the investigation of three-body correlations 
in a ID Bose gas. A magnetically trapped Bose-Einstein 
condensate of 87 Rb atoms was loaded into a deep 2D 
optical lattice, by which the condensate was divided into 
an array of independent ID systems. In the case of 87 Rb 
two-body losses are very small, and it was possible to 



model the decay in time of the number of trapped atoms 
in terms of one-body and three-body processes only. The 
rate of change of the total number of atoms N was written 

as 



dvK. 



3D: 



(2) 



where K\ and K^ are the rate coefficients for one- and 
three-body losses. By measuring the decrease of the num- 
ber of atoms as a function of time it was found that the 
rate coefficient in the ID gas was reduced consider- 
ably, by a factor of about 7, compared to its value for 
a 3D condensate. The measurements were carried out for 
a particular value of the dimensionless coupling constant 
7 (« 0.5) characterizing the strength of the ID correla- 
tions. Given that K 3 is proportional to the three-body 
correlation function 33, the observed reduction was inter- 
preted as an effect of the reduced dimensionality on g$. 

In the following we shall therefore focus on calculating 
53 for a ID Bose gas. Our treatment is based on the 
Lieb-Liniger (LL) model jj| for which the dimensionless 
coupling parameter 7 is given in Eq. (0 below. In the 
thermodynamic limit, where the particle number N and 
the size of the system L tend to infinity while the density 
n remains constant, 



= N/L — const, N,L — ► 00, 



(3) 



an expression for 33 has so far been obtained only for 
small and large values of 7 (3 ■ Here we report an explicit 
expression for #3, valid for all 7. Its derivation employs 
an integrable lattice regularization of the LL model jfl, |j| 
together with conformal field theory. The calculations 
are quite involved and will be presented elsewhere @|. 
The answer is given by Eq. I|17ll in terms of two thermo- 
dynamic parameters of the LL model: the second- and 
fourth-order moment, Eq. jlfijl . of the quasi- momentum 
distribution function. The quasi-momentum distribution 
function is the solution of the Lieb equation 11411 , which 
is a linear integral equation. 

Within the LL model the Hamiltonian for N identical 
bosons of mass m is taken to be 
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The interaction constant c > has the dimension of in- 
verse length. The dimensionless coupling strength 7 is 
given by 



7 = c/n. 



(5) 



The LL model allows for an exact determination of the 
eigenfunctions and eigenvalues of H. In the 7 — > 00 
limit of the model, known as the Tonks- Girardeau (TG) 
gas 0, the eigenvalues are the same as that of a free 
fermion gas. 

While the eigenfunctions and eigenvalues are known 
exactly for the LL model, its correlation functions in 
general and 53 in particular have been much less investi- 
gated. The reason for this is that the eigenfunctions are 
very complicated for a general value of N. We exhibit 
these below and calculate g-$ analytically for the simplest 
possible case N = 3. 

The eigenfunctions ^(xi, x%, . . . , xjy) of the Hamilto- 
nian iQJ must be symmetric under any interchange of 
coordinates. Since the particles interact via a delta- 
function potential, the derivatives of \& jumps when two 
particles approach each other, while when all dif- 
ferent, \t satisfies the free-particle Schrodinger equation. 
The explicit expression for VP was obtained in Ref. 0] and 
given in Ref. 10] in the following form 



N 



9(xx,x 2 ,...,x N ) = C^(-1) [P] exp<^ iJ2 k P 3 



x Yli k Pj ~ k Pi ~ «csgn(ajj - xi)], (6) 

3>l 

where P is a permutation of N numbers and [P] is the 
parity of the permutation. The possible values of the 
quasi-momenta kj are determined by the boundary con- 
ditions imposed on the system. For periodic boundary 
conditions (that is, for TV particles placed on a ring of 
circumference L) kj are solutions to the following system 
of coupled nonlinear equations, called the Bethe equa- 
tions, 



exp {ikjL} 



ki + ic 



,N. (7) 



The calculation of the normalization constant C in 
Eq. © is a nontrivial task. For the periodic bound- 
ary conditions a closed expression for C was suggested 
by Gaud in ( see fll] | and references therein) and proved 
in Ref. [lif : a detailed discussion can be found in 
Ref. 0| . One has now all the ingredients for calculat- 
ing gs = 33(7, N) for the ID system. Written in the 
first-quantized form, 93(7, N) is 



(JV-3)! 



dX\V(Q,0,0,XA,...,x N )\* 



(8) 
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Figure 1: The three-body correlation function (73, normalized 
to its value in the absence of interactions, as a function of 7. 
The dashed lines are for particle numbers JV = 3, 4, 5, and 6 
(from top to bottom), and the full line is for N = 00. 



where dX = dx^ . . . dx n and ^ is given by Eq. ©. Note 
that the average (...) in Eq. is taken over the ground 
state of the system; equivalently, it is the ground-state 
wave function which enters Eq. (JHJ. This is ensured by 
selecting the proper set of the quasi-momenta kj among 
all the solutions of the Bethe equations @, 

All quasi-momenta kj are equal to zero for the non- 
interacting system, 7 = 0, and the wave function ffjjj is 
uniform in space, \P = (1/VL) n , which implies that 



g 3 (0,N) = n 3 N(N -1)(N -2)/N 3 



(9) 



We now use the wave function Q to get #3(7, JV) for 
N = 3 analytically, and for N = 4, 5 and 6 numerically. 
The quasi-momenta in the ground state of the system 
with TV = 3 obey fci = ~k 3 and &2 = 0. It is convenient 
to write the Bethe equation for k\ in the form 

A = 2tt - 2 arctan(A/37) - 2 arctan(2A/37), (10) 

where we have introduced A = k\L and the princi- 
pal branch for arctangent has been chosen: — ir/2 < 
arctan(a;) < tt/2. The solution to this transcendental 
equation grows monotonically from A ~ 3^/7 for 7 <C 1 
to A ~ 2n — 471-/7 f° r 7 ^ 1- The resulting answer for 
ff3 (7,3) = 6|vE-(0,0,0)| 2 is 



53(7,3) = 
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3 7 3 L 3 [1 + 4(3A 2 + 2A + B + GAB)] ' 



(11) 



where A = 3 7 /(9 7 2 + A 2 ) and B = 3 7 /(9 7 2 + 4A 2 ). The 
ratio 53(7, 3)/g3(0, 3) is plotted in Fig. ^ It decreases 
monotonously with increasing 7, from 1 — 297/24 for 
7 < 1 to 512tt 6 /2437 6 for 7 > 1. In Fig. □ we also 
plot 33(7, N) calculated for N — 4, 5 and 6. Note that 
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the magnitude of the slope at the origin (7=0) increases 
with increasing number of particles, and when N — > oo 
it approaches infinity in agreement with the Bogoliubov 
result given in Ref. 



53(7, 00) /n 3 ~ 1 - 6^/7/ 7T, 7- 



0. 



(12) 



We have obtained an analytic expression for 33 in the 
Tonks-Girardeau limit, 7 — > 00, for arbitrary values of 
iV by generalizing the result Eq. (16) of Ref. [5j, valid in 
the thermodynamic limit, to the case of a finite number 
of particles. The result is 



53 (7, N) 16^ 6 (iV 2 -l) 2 (iV 2 -4) 



157 6 



N 6 



, 7^00. (13) 



We now present our main result: the exact expression 
for the three-body correlation function 33(7,00) in the 
thermodynamic limit of the Lieb-Liniger model J3J. 
The system of Bethe equations (0 reduces in this limit 
to the linear integral equation called the Lieb equation: 



a(k) 



1 



, 2aq(g) 
a 2 + (k — q) 2 



2vr j_ t 

where a is an implicit function of 7 : 
a = 7 J dk a(k). 
In terms of the moments e m defined by 
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m+l 



dkk m a{k), m = 2,4. 



(14) 



(15) 



(16) 



the exact expression for 33(7,00) is 



53(7,00) 



27 



f 4 
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1+ 1 
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3e 2 e 2 



where e' m is the derivative of e m with respect to 7. 

Equations l(T3jl - i(T7t form a closed set of equations de- 
termining 33(7,00). We plot the result in Figs. and EJ 
For a numerical solution of the integral equation ifHjl we 
have used the Mathematica package NISolve 0. InFig.Q] 
we compare the normalized function 33(7,00)733(0,00) 
with that calculated for N = 3, 4, 5, and 6. One can 
see that the convergence to the thermodynamic limit 
is quite slow. In Fig. we compare 33(7,oo)/n 3 with 
the experimental data from Ref. [!{]. The average value 
of 7, measured in [l|, 7m ~ 0.45, and the correspond- 
ing value of 33(71^,00) » 0.14 are shown as the dot in- 
side the box which represents the experimental uncer- 
tainty according to Ref. Q: 0.34 < 7 m < 0.65 and 
0.05 < 33(7 m ,oo) < 0.23. The dashed lines show the 
asymptotic expressions given in Eqs. 112t and 111 -ill for 




Figure 2: Double-logarithmic plot of the three-body corre- 
lation function (73(7,00) as a function of 7 (full line). The 
dashed lines indicate the asymptotic results Q given in 
Eqs. 11211 and Hl.'ifl for small and large 7, respectively. The 
measured value is shown as the dot inside a box indicating 
the experimental uncertainty Q. 



small and large 7, respectively. Evidently, the asymp- 
totic expressions do not account for the observed value, 
and the exact expression Q17JI is needed. 

Eq. itTTIl gives g3(-f,oo)/n 3 ~ 167r 6 /157 6 when 7 — > 
00, thus reproducing the asymptotic expression Eq. l|T3jl 
taken at TV = 00. To check this is a straightforward task 
since Eq. 114t admits a regular perturbative expansion 
with 1/a being a small parameter. The opposite limit, 
7 — > 0, of Eq. ijTTjl is much more difficult to analyze since 
the kernel in Eq. 11 411 becomes singular when a — > 0. 
Beyond the leading order, the results were obtained only 
recently [l^j : 



a(k) ~ v 7 ! - k 2 /2na + f(k), 







(18) 



where f(k) is written explicitly in Ref. [14j |. see Eq. (4.11). 
For the moments Eq. ltT6|l one gets e 2 — 7(1 — 4^/7/3^) 
and £4 ~ 27 2 (1 — 44^/7/15^). With these expressions for 
the moments, expansion ltl2|l is reproduced. 

A useful representation of our result l|17H is the follow- 
ing approximate form 



1 - 67r _1 75 + 1.26567 - 0.29597^ 

1 -0.22627-0.198172 
0.705 - 0.107 7 + 5.08 • 10~ 3 7 2 
1 + 3.417 + 0.903 7 2 + 0.4957 3 ' 
16tt 6 9.43 - 5.407 + 7 2 
15 7 6 89.32+ 10.197 + 7 2 ' 



0<7<1, (19) 
1<7<7, (20) 
7 < 7 < 30. (21) 



The relative error of this approximation does not exceed 
2-10 -3 for all values of 7 in the interval < 7 < 30, which 
we expect to be experimentally relevant. Note that the 
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asymptotic result 11 -til for N = oo is nearly a factor of 
two larger than the exact result for 7 = 30. 

The formal derivation of Eq. ifTTjl is quite lengthy and 
will be given elsewhere @|. Here we shall only sum- 
marize the main steps: We are interested in the ther- 
modynamic limit of the model, defined by Eq. Q). In 
this limit it is natural to use a field-theoretical approach 
rather than the first-quantized formalism we employed 
for the few-particle systems. The Lieb-Liniger model is 
completely integrable. In the field-theoretical language 
complete integrability implies the existence of an infinite 
family of conserved currents J n , n = 0, 1, 2, The in- 
tegrals over space of these conserved currents commute 
with each other, including the Hamiltonian <@J. There 
is thus an infinite set of mutually commuting operators 
J n , n = 0,1,2,..., for which the eigenfunctions and the 
spectrum are known exactly. A procedure for generating 
J n is discussed, for example, in Ref. The first three 
currents generated by this method are Jo = J\f, J\ = V, 
and J2 — H, where TV, V, and H are the number density 
of particles, the momentum density, and the Hamiltonian 
density, respectively. These three currents exist for non- 
integrable models as well, while the currents J n for n > 2 
are specific to the LL model. To illustrate our approach 
we introduce the classical expression for the current J± 
(its quantization will be discussed in the next paragraph): 

+ c(9 s #*) + 2c 2 (* 3 )** 3 . (22) 

The quantum-mechanical version of the last term in Eq. 
(12211 contains ^(x) 3 ^(x) 3 in which we are interested. 
When averaged over the ground state, this term gives 
the desired correlation function g 3 , Eq. Q. Some of 
the remaining terms in H22I1 can be eliminated by us- 
ing the Hellmann-Feynman theorem, which applies to all 
conserved quantities. This theorem was applied to the 
Hamiltonian of the LL model in Ref. to calculate the 
correlation function (w(x) 2 $(x) 2 ). In our case not all 
terms can be eliminated by use of the Hellmann-Feynman 
theorem. We have succeeded in obtaining additional 
identities || by employing the long-distance asymptotics 
of the correlation functions of conserved currents known 
from the conformal limit of the theory. 

The quantization of the classical current Eq. i122li 
is far from being straightforward. A naive quantum- 
mechanical version of the expression (1*221) contains terms 
involving derivatives of field operators and . Such 
terms are not well-defined in general. An explicit demon- 
stration of this was considered in the paper [1J|, where 
it was shown that the expectation value (iff* (0)<9 3 *I' (x)) 
is discontinuous at x = 0. The reason for this behav- 
ior can be seen from the structure of the Bethe-ansatz 
wave function ©, which has a cusp whenever the co- 
ordinates of any two particles coincide. This does not 



cause difficulties when one considers the momentum op- 
erator and the Hamiltonian itself, that is, when one works 
with the operators * f (x)d x $(x) and * t (x)5^(x). How- 
ever, for other operators we encounter difficulties which 
we have not been able to resolve within the continuum 
model l@J. Therefore we have used an integrable lattice 
model, the so-called g-boson hopping model 0,0; which 
reduces to the Lieb-Liniger model (@J in the continuum 
limit. Within this model all higher (quantum) integrals 
of motion are defined unambiguously. 

In conclusion, we have obtained an exact expression 
for the three-body correlation function Eq. within the 
Lieb-Liniger model of a one-dimensional Bose gas. The 
explicit expression is given by Eq. ifTTIl and Eqs. iflfljl - 
ll2l)l . and compared in Fig. [2] with a recent experiment. 
The finite-size effects are explored and a comparison with 
the thermodynamic limit is given in Fig. We expect 
that our results will stimulate further experimental in- 
vestigations of the local correlations in low-dimensional 
trapped atomic gases. 
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